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Abstract 

The onset of condensation of hard spheres in a gravitational field is 
studied using density functional theory. In particular, we find that the 
local density approximation yields results identical to those obtained 
previously using the kinetic theory [Physica A 271, 192, (1999)], and a 
weighted density functional theory gives qualitatively similar results, 
namely, that the temperature at which condensation begins at the 
bottom scales linearly with weight, diameter, and number of layers of 
particles. 

1 Introduction 

In a recent paper, one of the authors (DCH) [1] proposed that hard spheres 
in the presence of gravitational field, undergo a condensation transition, 
and identified the transition temperature, Tc, as a function of external pa- 
rameters, i.e.: 

Tc = mgDfi/fXo (1) 

where m and D are the mass and diameter of the hard spheres, /i is the 
dimensionless layer thickness at T = 0, and Hq is a constant that reflects the 
particular manner in which a system packs upon condensing. It was noted 
that there exists a critical temperature Tc below which the total number 
of particles is not conserved. This is the temperature at which the density 
at the bottom layer becomes the close-packed density. Now, since the hard 
spheres cannot be compressed indefinitely, if the temperature is lowered be- 
low Tc, then the first layer should remain at the close-packed state, while 
the particles at the second layer try to compact themselves and thus crys- 
tallize. The crystallization then proceeds upward from the bottom layer as 
the temperature is lowered. This picture was later confirmed by Molecular 
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Dynamics simulations [2] for mono-disperse hard spheres and was extended 
to the segregation of binary mixtures of hard spheres of different mass and 
diameters [3]. 

In the original work [1], the Enskog kinetic equation was used to obtain 
the density profile of hard spheres under gravity. However, in an attempt to 
solve a highly nonlinear integro-differential kinetic equation, it was assumed 
that the equilibrium velocity distribution function, /(r, v), factorizes into 
a product of space and velocity dependent parts, i.e. /(r, v) = G'(r)0(v) 
and further that the functional form of 0(v) is Gaussian, which should be 
valid for elastic hard spheres. The factorization assumption is an equilib- 
rium ansatz, which states that the configurational statistics are separated 
out from the kinetics when the system is at equilibrium, so that all the equi- 
librium quantities can be obtained from the configurational integral of the 
partition function. The equilibrium state is then the configuration that min- 
imizes the free energy. Therefore, we find it necessary to obtain the results 
of Ref. [1] by the variational method. Indeed, we will show in this paper 
that essentially identical results follow from an application of density func- 
tional theory (DFT) [5] for liquids to the problem. We will first employ the 
simplest form of the density functional theory known as the local density 
approximation (LDA), which assumes that the range of inter-particle inter- 
action is much smaller than the typical length scale on which p(r) varies [6]. 
We will show that the LDA and the Enskog theories are in fact identical, 
so that in both methods, the condensation temperature is defined as the 
temperature at which the sum rule breaks down. Next, we will analyze the 
problem with a simple weighted density approximation (WDA) [6-9] , which 
takes into account the local variation of the density function. In this ap- 
proximation, the condensation temperature is defined as the temperature at 
which the volume density at the bottom layer reaches the maximum allowed 
value. In this approximation, microscopic information is preserved in the 
density profile; notably the formation of a crystal shows up in the density 
profile as oscillations. The peak to peak distance of this oscillation is ap- 
proximately the particle diameter. We will demonstrate that the results of 
both analyses present a picture identical to those presented in Ref. [1]; in 
particular we will show how the value /iq that appears in Eq. (1) depends on 
the approximation. 
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2 Local Density Approximation 



The essence of the LDA is to assume that the system may be divided into 
small pieces of nearly constant density and then to treat each piece as though 
it were part of a homogeneous system [6] . Under these assumptions one may 
write a free energy functional: 

Flda[p\ = J (irp(r)V^(p(r)) + J dvp{v)UU^), (2) 

where '^{p{r)) is the Helmholtz free energy per particle in the absence of 
an external field and C/e^t is the potential energy per particle due to an 
external field such as gravity. Minimization of this functional under the 
global constraint that the number of particles is given by 

I drp{r) (3) 

should yield the desired density profile. To be more specific, we define 
variables for hard spheres of mass m confined in a d dimensional volume 
V = L'^~^H with L'^~^ being the cross sectional area of a {d — 1) dimensional 
plane and H being the height of the container along which the gravitational 
field is acting. The Helmholtz free energy per particle consists of two terms, 

^{p) = ^id{p) + ^exc{p), (4) 

where ipid is the ideal gas contribution, 

V^,,(p) = T(log(AV) - 1). (5) 

Note that ipid = —T\og{z^)/N with the single particle partition function 
z = V{27rmTY'/N\ [10] and that we have redefined the thermal wavelength 
A = (27rmT)^/^. Next, ipexc is the excess contribution to the free energy 
is due to the configurational integral coming from the interactions among 
particles and is in general written as the integral: 




The above equation can be derived from the thermodynamic relation, P = 
-(If )t with the chain rule: (^)r = ^(|-)t. Note that P/pT - 1 is the 
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virial sum. Since gravity acts along the vertical direction U^xt = rngz^ 
and the transverse degrees of freedom can be integrated out to yield the free 
energy functional per unit area: 



Flda\p] _ 
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oo 



_ roo 

= F\p] = / dz p{z)^id{p{z)) + 
Jo 

dz p{z)lpexc{p) + 
roo 

mg / dzp(z)z, (7) 
Jo 

where A = is the cross sectional area in the x — y plane. Minimization 
of the functional under the constraint Eq. (3) yields an equation for the 
density profile p: 

^ = T\og(A'p(z)) + ^P^M + -rmgz^X (8) 

where we have introduced a Lagrange multiplier, A. Defining = pD^ and 
C,= z/D, A should be determined by the sum rule: 

r(t>o roo 

//0e= / d<t>a<t>)- / dcm. (9) 

Jo Jo 

where 0c is the close-packed density, p is the number of layers of particles 
in the system at T = 0, and 00 is the density at the bottom layer. Note 
that the particular shape of the density profile will depend on the functional 
form of the pressure P or, equivalently, on the functional form of the excess 
free energy, if^exc- One may use the Enskog pressure for hard disks or a hard 
sphere equation of state given by a functional form: 

P = pT\l + ^pD\{p)l (10) 

where is the pair correlation function evaluated at contact {r — D), and 
where 7 = f when d = 2 and 7 = ^ when d = 3. Then 

,^djy 
p' 

Substituting this form of ipexc into Eq. (8) and taking the derivative with 
respect to z generates, in our non-dimensional variables, the differential equa- 
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which is precisely the result obtained in Ref. [1]. Thus the equivalence 
between the LDA and Enskog theory has been shown, and the constant /iq 
that appears in Eq. (1) can also be derived by the density functional theory 
in the local density approximation. 

To conclude this section we cite some results of the LDA/Enskog theory. 
Note that for the hquid phase, the density profile 0(C) in Eq. (9) is a mono- 
tonically decreasing function of the height ( with its maximum vahie at the 
bottom. Further, the maximum density 0o is a function of temperature, too, 
with the upper bound 0o < <Pc- So, the right hand side of Eq. (9) can be 
written as f{4>o)/P, where /3 = mgD/T. The particular form of the function 
/(0o) depends on the approximation. Ref. [1] (Eq. (15c) and (16c)) gives the 
functional forms of 0o in 2d using the Ree and Hoover correlation function 

[11] 

(1 - q;i0 + q;20^) , . 

^^'^^ ^ (1 - a<i>r ' ^''^ 

a = 0.489351 7r/2, 
ai = 0.196703 7r/2, 
a2 = 0.006519 7rV4. 

and in 3d using the Carnahan-Starling equation of state: 

P _ {I + T] + rj^ - Tif) 

(1-77)3 



(14) 



where we have defined the volume fraction rj — ^D^p = |0 in 3d and as 
rj = \D^p = |0 in 2d. For completeness, we reproduce these. First, the 2d 
result: 



1 ,2 , C3(/), 



fMuH = (1 + ^2)00 + ^Ci0^ + ^Y^^ (15) 
Ci ( 1 ^\ _^ C40O 



a \ (1 — achn) / (1 — ac 



\2 



with ci = TTcia/a^ ~ 0.0855, C2 = -{n /2){ai/ - 2a2lo?) ^ -0.710, cg = 
— C2, and C4 = (7r/2)(l/Q! — axjo? + aijo?) ~ 1.278. Next, the 3d result: 

(1 - Q;0oj (1 - Q;0oJ 
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where in this expression a — n/Q. Note that in both 2 and 3 dimensions, 
/(0o) is a monotonically increasing function of 0o. Hence, it has a maximum 
at 00 = 0c, the value of the close-packed density. Since j3 or equivalently T 
and the layer thickness n are arbitrary control parameters, the sum rule, Eq. 
(9), breaks down when T <Tc, where 

A*0c = fmaxTc/mgD = iioTjmgD (17) 

While the scaling of the critical temperature displayed in Eq. (1) is 
independent of the particular equation of state used in the calculation, the 
maximum value of /(0o), A^o, depends on the functional form of the density 
profile, or equivalently, the pressure. Using the two approximations above, 
and taking the maximum densities as r]c = f \/2 ~ 0.74 in 3d, and r]c = 
tt/{2V3) ^ 0.91 in 2d, we find 

f,ORH = 111-31 (2rf) 

l^ocs = 152.34 {3d). (18) 

At the level of Enskog approximation, no is quite sensitive to the density at 

the bottom, 0o. At this point, we find it appropriate to mention the point 
made by Levin [16] that reliable information about the fiuid-solid coexistence 
cannot be obtained by the LDA, because of its inability to include the density 
variations in a highly structured phase (solid) . When the Enskog approxima- 
tion breaks down, one has to either abandon the approximation and search 
for a better one, or modify the approximation by removing the unphysical re- 
sults. In the original paper [1] , the latter approach was taken, namely based 
on physical grounds, the crystal regime was replaced by a constant average 
density, a Fermi rectangle, and the fiuid regime was then fit to the Enskog 
profile, which is linked to the Fermi rectangle at the hquid-sohd interface. 
While the proportionality constant fio in Eq. (1) obtained this way seems to 
overestimate, and thus while the Enskog equation fails to locate the precise 
point of the liquid-solid transition, the prediction of its existence and the 
scaling relation between the critical temperature and external parameters 
(Eq. (1)) seem to remain true. More elaborate approximations that do take 
into account the local variations in the structured phase yield substantially 
lower values for /xq (see Fig. 2), which are somewhat close to the values ob- 
tained by a Mean Field theory [2]. In order to show the dependence of fiQ on 
approximation, we also compute in it 3d by the Percus-Yevick compressibility 
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form of the equation of state: 

P 1 + 77 + 77^ 



pT (1 - 77)3 
which yields equally high values for /iq, 



pT (1 - 77)^ 



yields 



(19) 



(1 - acpo) (1 - a0o) 

00 

l^opYc = 185.19. 
The slightly different form, namely the virial form 

P 1 + 2// + 3if 



(21) 



(l-a0o) (l-"0o)^ 

A^opyy = 86.63. 

We further point out that the breakdown of the sum rule is due to the 
fact that the pressure has a singularity at 77 = 1, and thus it has a finite value 
at the close-packed density 77c, which is necessarily less than one [10]. If one 
uses the lattice gas pressure [4], 

P = -Tlog(l-p), (23) 

which has a singularity at p = 1. then the condensation temperature is zero, 
and the density profile is given by the Fermi function [5]: 

p{z) = 1/(1 + exp(7n^(^ - p)/T)) (24) 



3 Weighted Density Approximation 

The essence of the WDA, as introduced by Tarazona [6,7] and Curtin and 
Ashcroft [8] is to recast Eq. (2) , the general form of the free energy functional. 
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as 

Fwda\p] = j dvp{v)%pid{p{r)) + j dvp{v)ip^^^\p^{Y)] 

+ j drp(r)Ue,t(r), (25) 

where ipexc[PwiT^)] is now a functional of p(r), depending on p(r) through the 
weighted average of the density given by 

p^(r) = Jdh'w{\v-v'\)p{r'), (26) 

where w(|r — r'|) is an appropriately chosen weighting function. Following 
Tarazona [6], we choose 

where © is the unit step function, i.e., we replace the local density p(r) with 
its average over a sphere of radius equal to the particle diameter D. Because 
we assume planar symmetry, i.e., independence in the x and y directions, we 
may integrate out the transverse degrees of freedom and write explicitly the 
integral above as a one dimensional integral ior z > D: 

One needs to be careful near the bottom layer z = 0, namely, for < z < D. 
In this case, the weighting cannot be done over a sphere a radius D, because 
of the infinite potential at 2; = 0. We propose to carry out the weighting 
over that part of the sphere of radius D and centered at z above the z = 
plane. Thus, the normalization factor, C, that is, the volume over which the 
integration is performed, is no longer C = j_^{TT{D'^ - z'^)dz' = A-kD^/^, but 
instead is C = J!^,{tt{D^ - z'^)dz' = TrffD^ + D^z - ^z^]. Hence, for z < D, 

poo 

Pw{z) - [|g3+fl2,_i,3i dz'p{z'){D^-{z-z'f)x 

e{D-\z-z'\). (29) 

As before, we need to extremize the free energy functional under the 
global constraint on particle number, so we again use the method of Lagrange 
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multipliers and functional differentiation. Performing the minimization of the 
free energy functional (Eq. (7) with p in the excess term replaced by (28)), 
we find the following equation must hold: 



POO 

Tln(AV) + i^exc{pM)+ dz'p{z') 

JO 



,Jlpexc{pw{z')) 



5p{z) 

+mgz + A = 0. (30) 



We write explicitly the integral term in the equation above: 

.xci^Pw ( 



^B{z,z')T{z'), (31) 



iiz' > D and 



^ " dpM') ' ^^'^ 

B{z, z') = [d^ -{z- z'f) QiD -\z-z'\), (33) 
= [|2^3 + _ 1^.3] (35) 

iiO<z' <D. 

The integral equation for p{z), Eq. (30), is highly nonlinear and complex. 
Thus, it requires numerical solution. We choose to solve Eq. (30) using the 
Carnahan-Starling equation of state, Eq. (14), so that 

A{z') = 



p^{z') (l-^D^p^z')) 

' (36) 



p^{z') [l-^D'^p^{z' 

For a given choice of A we iterate Eq. (30) until the iteration converges to a 
unique profile. The integral of the profile (Eq. (9)) determines p, effectively 
the depth of the unexcited sample, so for fixed m, g, D, and T, we tune A to 
control the number of particles. 
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Figure 1: (a) The volume density 77 as a function of the dimensionless height 
z/D at different temperatures as calculated from numerical solution of Eq. 
(30) for given set of m, g, D, and fi. The tail extends as the temperature 
increases. At high temperature T, the density profile is a monotonically 
decreasing function of z/D, but as the T approaches the critical temperature 
Tc, oscillations develop near the bottom, indicating layer formation, (b) The 
dotted line is a magnified view of the topmost curve in Fig. la, while the 
solid line is the prediction made by the LDA/Enskog theory for the same 
system. 



We find that at high temperatures, the profiles obtained using the WDA 
match very well the profiles obtained for the same set of parameters using 
the LDA/Enskog approach. But as we lower the temperature of the system, 
particles at the bottom begin to compact themselves, and the crystaUization 
sets in. One of the notable features of this weighted density functional ap- 
proach is that formation of the crystal can be captured in the density profile, 
in particular oscillations in the density profile appear near z = 0. Fig. la 
shows the development of these density peaks for a representative system at 
three different temperatures above Tp. The peak to peak distance of this 
oscillation is slightly greater than the diameter of the hard sphere. Fig. lb 
is a closer view of the density profile for the coolest of these systems (dot- 
ted line). This figure also plots the LDA/Enskog result for the same system 
(solid line); it agrees well with the WDA profile for large z/D, but cannot 
reflect the rapid oscillations in density which occur near the bottom of the 
sample. 

With sufficiently low temperature, the bottom-most peaks in the profile 
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Figure 2: The dimensionless condensation temperature, Tc = TjmgD, 
is plotted against the dimensionless layer thickness /i. The slope is l/no. 
Squares are the result for d = 2, circles for d = 3. 

approach (and even exceed) the physical hmit of = 1. Note that the maxi- 
mum 1] for close-packing in 3d is 0.74. However, in our numeric solution, we 
have chosen the lattice spacing, Az = D/256, with D the particle diameter. 
This choice is to guarantee precision in the solution of the integral equation 
Eq. (30). Hence, even though the physically relevant hmit for the volume 
fraction r] in 3d is 0.74, in our numerical solutions, 7] must approach one at the 
close-packed density. With this modification, we define the temperature at 
which f] first reaches this physical limit at 2; = as the critical temperature, 
Tc. In Fig. 2, for a given set of rrijgjD, we have plotted the dimensionless 
critical temperature Tc = Tc/mgD as a function of the initial layer thickness, 
ji. The numerically determined value from the slope for the constant yUo is 
/io = 6.10 in 3d. We have performed an analogous WDA calculation in 2d 
using the Ree and Hoover correlation function x(0), Eq. (13). The data from 
this calculation also appear in Fig. 2. They yield //q = 3.52 in 2d. Both the 
2d and 3d WDA results are smaller than those obtained by the LDA/Enskog 
approach, and we discuss this next. 

As we have discussed, in the LDA/Enskog approach, the value of /Uq de- 
pends on 00, the density at the bottom, and is identical to the function /(0o)- 
In all the approximations we have used in this work, /(0o) (see Eqs. (14, 15, 
20, 22)) is a function very sensitive to 0o for 0o near close-packed values, i.e. 
for (j)Q>l. Fig. 3a illustrates this sensitive dependence. Molecular dynam- 
ics simulations in two dimensions [12] have shown that 0o at T < Tc varies 
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Figure 3: (a) The value /Xq = /(0o) as a function of density at the bottom 
of the sample, as calculated in the LDA/Enskog theory. The solid curve 
is for 2d using the Ree and Hoover value of x(0), Eq. (15). The remaining 
curves are for 3d: dotted, Percus-Yevick compressibility form, Eq. (20); 
dashed, Carnahan-Starling, Eq. (16); long dashed, Percus-Yevick virial form, 
Eq. (22). The arrows on the 2d curve indicate the range of fxo calculated 
with the LDA/Enskog theory from molecular dynamics simulation values of 
00- (b) An expanded region of Fig. 3a, showing only the 2d Ree and Hoover 
curve (solid) and the 3d Carnahan-Starling curve (dashed). See the text for 
further discussion. 

widely. For one set of simulations using 10^ hard disks with = 20, defects 
in packing lead to 0o occupying the range 1.00 < 0o < 1-14, with higher den- 
sities occuring at lower temperatures. (Note that for square packing in 2d, 
00 = 1, while for triangular packing, 0o = 2/\/3 ~ 1.155.) In LDA/Enskog 
theory, this range of 0o leads to 21.76 < Hq < 90.33 (the arrows in Fig. 3a 
indicate this range), the large range due to the sensitivity of f{(po), while the 
WDA theory presented here gives a smaller value, //q = 3.52 for the 2d cal- 
culation using the same equation of state. Though the discrepancy between 
the WDA and LDA/Enskog results seems large, the source of the discrep- 
ancy is easy to identify. To do so, we first remark that any density profile 
derived through the WDA for a given system (m, g, D, jji, and T) at tem- 
perature below Tc (WDA) differs from the LDA/Enskog profile for the same 
system appreciably only near the bottom of the sample, where the WDA 
profile exhibits oscillations and the LDA/Enskog profile is a monotonically 
decreasing function bounded between the peaks and troughs of the WDA 
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profile (see Fig. lb). At the temperature when the bottom-most density 
peak in the LDA profile reaches its maximum value (f>max = ^Vmax = ^, i-e. 
at Tc (WDA), the LDA/Enskog profile for the same system will have a much 
smaller maximum. In our work in 2d, at Tc (WDA), the LDA/Enskog profile 
has 00 = -Vo — - 0.55 0.70, below the square packing value, while in 3d, 
at Tc (WDA), the LDA/Enskog profile has 0o = f r^o = f 0.42 ^ 0.80, below 
the simple cubic packing value. If we use these values to compute fiQ in the 
LDA/Enskog approach as in Fig. 3b, we get fiQ = 3.56 in 2d and Hq = 6.20, 
consistent with our determination of /xq from the slopes of the lines in Fig. 2. 
We see then that Tc (WDA) is in general higher than Tc (LDA/Enskog) for 
the same systems, and this is reflected in the lower value of /Iq in the former 
approach. 

Finally, we turn our attention to the question of whether the condensation 
phenomenon we are considering is a phase transition in the thermodynamic 
sense, i.e., whether condensation corresponds to a discontinuity in the flrst or 
higher derivatives of the free energy with respect to temperature. We address 
this question by focusing on the gravitational potential energy contribution 
to the free energy, Ug = mg /q°° zp{z)dz, which is proportional the center of 
mass < z >= zp{z)dz/ p{z)dz. First we show that in the LDA/Enskog 
theory, which is extended to temperatures below Tc by the assumption that 
the density in the frozen layers is given by (j) — (pc and that the density 
above the frozen layers is given by a vertically shifted LDA/Enskog profile 
[1], a kink in the center of mass develops at T — Tc, suggesting a first order 
transition. 

To do this we recall that the density profile 0(C) is given by the functional 
form: 




- /(0o) 

= mgD/T. Above T^ 



(37) 



< CiT) > 




(38) 



where 
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Now, we note that for T near Tc- 

4>o{T)^ct),-a{T-T,) (39) 

where a > 0. Then for any integrand G{4>) we can make the following 
approximation: 

r G{4>)d4> ^ f G{4>)d4> - a{T - T,)G{4>,). (40) 

Applying this to the above expression to the integral Ii and I2, we find that 
< C(^) > is linear in T with a quadratic correction. 

Below Tc, the density profile develops a kink ai — L. For ( < L, 
0(C) — 0c the close-packed density, and for C > L, the profile is given by the 
LDA/Enskog profile Eq. (37), and the thickness of the frozen layer is given 

by[i] 

L = M1-T/T,). (41) 
We now compute the center of mass < C{T) >■ 



(42) 



where 



I 



poo poo 

/ cmdc+L <f>{c)dc=h+Lh 

Jo Jo 



Hence, 



where 



I2 = Ml^-L). (43) 



0e//'^(l-^) + J 
c c 



Jo J4>c dcp 

and where 

A = j^\f{<t>c)-f{m'^d(t>. 
14 



It therefore follows that 



< C(T) >= I + AiT^ (44) 

where 

Ai = [-(^)^(^-f)]. 

The center of mass scales with temperature quadratically below T,. but lin- 
early just above Tc; thus, there is a kink in the center of mass and in the 
gravitational potential energy contribution to the free energy, giving rise to 
a first order transition. 

The scaling of < > with below Tc survives a modification of how we 
represent the frozen region. Suppose, the density in the frozen region is not 
represented by a uniform (p^ but is instead given by 

m-Y.pAc-Q (45) 

i 

where Q is the position of the center of hard spheres and pi is its peak density 
in the i-th row forming a crystal. This is a crude way to approximate the 
oscillations in the density profile due to the crystallization. Then, /i in Eq. 
(43) is replaced by: 

h= ^\mdC = Y.C^P^ (46) 

Jo . 

If Pi — 0c for all i, then, 

h = 0eEO = 0c[l/2 + 3/2 + 5/2 + .... + (2L-l)/2] 

i 

= 00^72 (47) 

which is the same result as that obtained by assuming the density profile is 
approximated by a Fermi rectangle. 

The WDA approach to the problem also yields results suggestive of the 
existence of a first order phase transition. It may be an artifact of the method 
that in the WDA solutions, as T is lowered below Tc as we have defined it, 
large density peaks whose maxima exceed the physical limit r] — 1 appear. 
Thus the WDA does not capture the exact distribution of material in the 
system. However, it is nonetheless suggestive to examine the dependence of 
the center of mass < C(2^) > oii T. Fig. 4 shows results for a representative 
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Figure 4: WDA calculation of center of mass < C > vs. T/mgD for a system 
with II — 10. The arrow indicates Tc and points to what may be a kink in 
the function. 

system in 3d with fi = 10 and whose critical temperature was found to be on 
the range 1.4 mgD < Tc < 1.5 mgD. An elbow, possibly a kink, is apparent 
in the vicinity of Tc, marking the onset of near linear behavior for T > Tc. 
We do not assert that this is evidence of a phase transition; we display this 
data merely to suggest that the existence of such a phase transition in the 
WDA approach is not inconsistent with our data. A different form for the 
weight function in Eq. (26) might yield a better result regarding the nature 
of the phase transition. 

4 Conclusions 

We have shown that the conclusion of the original paper [1], namely, that the 
scaling of the critical temperature at which hard spheres under gravity begin 
to form a solid is linear with their weight, their diameter, and the depth of the 
sample, necessarily follows from the simplest density functional theory for the 
problem (the LDA) and survives a richer density functional treatment using 
a WDA. Prudence requires us to note that our WDA for this problem did 
not include any sophisticated attempt to represent the crystal-fluid interface, 
something other researchers [13-15] working on similar problems have done. 
Doing so should likely give a more accurate quantitative picture than that 
presented here. 
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